classdef sfTuning <tuning
    
    properties
        units = 'cpd';
        calcMethod
        tuningParams
    end
    properties(Transient = true)
        sampleTuning
        sampleRange
    end
    methods
        function s = sfTuning(calcMethod,tuningParams,type)
            if ~exist('type','var')||isempty(type)
                type = 'sfTuning';
            end
            s = s@tuning(type);
            if ischar(calcMethod)
                ok = false;
                switch calcMethod
                    case 'rfModel'
                        % param needs KC,KS,RC,RS
                        if all(ismember(fieldnames(tuningParams),{'KC','KS','RC','RS'}))
                            ok = true;
                        end
                    case 'gaussian'
                        if all(ismember(fieldnames(tuningParams),{'optSF','sigmaSF'}))
                            ok = true;
                        end
                end
                if ~ok
                    error('sfTuning:inappropriateInput','calcMethod is wrong or params do not match calcMethod');
                end
                s.calcMethod = calcMethod;
                s.tuningParams = tuningParams;
            else
                error('sfTuning:wrongType','type needs to be a char');
            end
        end
        
        function out = sample(s,stim)
            switch s.calcMethod
                case 'rfModel'
                    FS = getSpatFreq(stim);
                    KC = s.tuningParams.KC;
                    KS = s.tuningParams.KS;
                    RC = s.tuningParams.RC;
                    RS = s.tuningParams.RS;
                    % needs to be scaled according to the maximum F1
                    Fss = logspace(log10(1/(10*RS)),log10(2/RC),100);
                    F1s = abs((KC*pi*RC^2*exp(-(pi*RC*Fss).^2))-(KS*pi*RS^2*exp(-(pi*RS*Fss).^2)));
                    factor = 1/max(F1s);
                    out = factor*abs((KC*pi*RC^2*exp(-(pi*RC*FS).^2))-(KS*pi*RS^2*exp(-(pi*RS*FS).^2)));
                case 'gaussian'
                    FS = getSpatFreq(stim);
                    out = exp(-(FS-s.tuningParams.optSF)^2/(s.tuningParams.sigmaSF^2));
            end
        end
    end
end